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Preface 


The  study  reported  herein  is  the  culmination  of  a  desire  which  X 
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adverse  pressure  gradients  in  several  previous  efforts  has  made  this 
study  very  self  satisfying.  When  faced  with  the  task  of  finding  a  thesis 
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my  understanding  of  current  numerical  techniques  and  the  state  of  their 
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interaction  of  a  shock  wave  with  a  turbulent  boundary  layer.  I  was  both 
surprised  and  pleased  to  find  that  Capt  John  Shea  had  accomplished  and 
reported  in  the  last  year  almost  exactly  what  I  wanted  to  undertake. 

The  effort  reported  herein  has  satisfied  both  of  my  personal  objectives 
I  by  extending  Shea's  numerical  study  of  normal-shock /turbulent  boundary 

layer  interactions. 
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Abstract 


'The  hybrid  finite  difference  code  developed  by  MacCormack  was 
applied  to  the  investigation  of  transonic  normal-shock  turbulent  boundary 
layer  interactions.  The  computations  were  performed  for  the  half  plane 
of  a  symmetric  two  dimensional  duct  by  establishing  a  symmetry  boundary 
condition  at  the  upper  boundary.  Both  first  and  second  order  center  line 
boundary  conditions  were  imposed  with  no  measurable  difference  observed. 

A  two-point  linear  extrapolation  of  the  primative  variable  was  unsuccess¬ 
fully  attempted  at  the  subsonic  outflow  boundary,  but  a  simple  zero  gra¬ 
dient  condition  gave  satisfactory  results  at  four  different  outflow 
boundary  positions  relative  to  the  shock  wave.  Numerical  results 

(M=  1.51,  1.40  and  1.30  Re  =  3-X-4Q  per  ft)  were  compared  with  the 

”C  a 

experimental  data  reported  by  Abbiss  and  East.  Even  though  the  data 
exhibit  three-dimensional  effects,  the  two-dimensional  computations 
show  agreement  within  approximately  10%.  The  differences  observed  in 
the  numerical-experimental  comparisons  were  all  consistent  with  expected 
three-dimensional  trends.  Although  not  conclusive,  the  potential  of 
adding  simple  three-dimensional  corrections  to  the  two-dimensional  code 
shows  promise  for  improving  the  experimental-numerical  agreement. 
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Nomenclature 


English  Symbols 


=i 

H,  h 

k 

M 

P 

R 

Re 

RL 

T 

tch 

u 


x,y 

Greek  Symbols 
a 

A 
Y 


Speec  of  sound 

Local  coefficient  of  friction 

2 

Specific  heat  at  constant  pressure,  for  air  6006  ft  /sec  R 

2 

Specific  heat  at  constant  volume,  for  air  4290  ft  /sec  R 
Total  energy  per  unit  volume 

2  2 

Specific  internal  energy  e  =  0.5  (u  +  v  )  +  e^/p 

2  2 

Enthalpies,  H  =  h  +  u  +  y 

2 

Coefficient  of  heat  conductivity 

Mach  number 

Pressure 

2 

Gas  constant,  1716  ft  /sec  R  for  air  R  =  c  -c 

p  v 

Reynolds  number 
Reference  length 
Temperature 

Characteristic  time,  (XT  ,  X  )/U 

Velocity  component  in  X  direction 
Velocity  component  in  Y  direction 

Reference  coordinate  system  with  x  parallel  to  wall 
(direction  of  mean  flow)  and  y  perpendicular  to  the  wall 


Local  flow  direction  measured  counterclockwise  from  the 
x  axis 

Finite  difference  differential  element 
Ratio  of  specific  heats,  Cp/cv 


ix 
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Boundary  layer  thickness 
Boundary  layer  displacement  thickness 
Viscosity  coefficient,  X  =  -2/3  y 
Molecular  viscosity 
Kinematic  viscosity,  y /p 


Density 


Shear  stress 


Subscripts  and  Superscripts 


Free  stream  or  unperturbed  condition 

Index  for  computational  mesh  in  x  direction,  i  =  1  at 
inflow  boundary 

Index  for  computational  mesh  in  y  direction,  j  =  1.5  at  wall 
Time  indices 

Total  or  stagnation  condition 
Derivative  with  respect  to  x,  y 


Computer  Code  Variables 


Column  index  for  parameter  storage  (I  fixed  for  any  x 
location) 

Outflow  boundary  value  of  I  index 


Row  index  for  parameter  storage  (J  fixed  for  any  y  location) 


Center  line  boundary  value  of  J  index 
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Introduction 


Background 

Green  (Ref  9)  and  Schofield  (Ref  23)  have  pointed  out  that  shock- 
boundary  layer  interactions  are  the  predominant  phenomena  in  many  aero¬ 
dynamic  problems  of  interest.  More  specifically,  normal-shock  turbulent 
boundary  layer  interactions  are  found;  (1)  in  transonic  and  supersonic 
inlets,  (2)  on  transonic  wings,  (3)  in  compressor  cascades,  (4)  on  rotor 
tips,  (5)  in  nozzles,  and  (6)  on  flow  over  control  surfaces  and  high  lift 
systems  (see  Figure  1) .  Although  they  have  been  a  primary  subject  of 
study  for  over  thirty  years,  very  little  progress  had  been  realized  in 
the  prediction  of  strong  interactions  until  the  last  three  to  five  years. 
The  lack  of  progress  was  due  to  the  complex  nature  of  the  flow  which  is 
produced,  and  the  nonlinearity  of  the  differential  equations  required  to 
model  the  flow.  The  application  of  todays  high  speed  digital  computers 
and  modern  numerical  techniques  (Refs  2,  5,  6,  14,  16,  20,  21  25-31)  shows 
great  promise  of  providing  useful  engineering  solutions  to  many  interaction 
problems . 

First  a  brief  discussion  of  the  physical  phenomenon  is  provided. 

The  viscous  flow  over  any  solid  surface  results  in  a  boundary  layer  near 
the  surface.  If  there  were  no  boundary  layer,  a  shock  could  theoretically 
be  generated  at  or  reflected  from  the  solid  surface.  However,  in  the 
presence  of  a  boundary  layer,  which  is  always  physically  present,  an 
abrupt  pressure  rise,  such  as  that  caused  by  an  impinging  shock  wave, 
cannot  occur  at  the  solid  surface.  The  subsonic  part  of  the  boundary  layer 
is  unable  to  support  such  a  sharp  pressure  change.  Consequently,  the  pres¬ 
sure  rise  is  spread  upstream  in  the  subsonic  region  of  the  boundary  layer. 
This  pressure  rise  causes  a  sudden  increase  in  the  boundary  layer 
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Figure  1 
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Sources  of  Normal— SYiock/Turbulent  Boundary  Layer  Interactions 
in  Nature 


displacement  thickness  and  precipitates  compression  waves  in  the  super¬ 
sonic  region  of  the  boundary  layer.  The  pressure  rise  is  therefore 
spread  smoothly  over  a  distance  on  the  wall  even  though  it  is  more  sharply 
defined  across  the  shock  wave  in  the  outer  flow.  The  distance  over  which 
the  pressure  rise  is  spread  is  on  the  order  of  the  undisturbed  boundary 
layer  thickness  unless  separation  is  induced.  For  cases  where  the  boun¬ 
dary  layer  separates  the  affected  region  can  be  much  larger  (Ref  5:4-5). 
This  shock  wave  configuration  is  commonly  called  a  lambda  shock,-  and  the 
point  where  the  shock  and  the  compression  waves  meet  is  known  as  the 
bifurcation  point. 

Shock  boundary  layer  interactions  are  usually  divided  into  two  major 
classes  commonly  denoted  as  strong  and  weak  interactions.  The  classifi¬ 
cation  is  dependent  on  whether  or  not  boundary  layer  separation  effects 
on  the  global  flow  result  from  the  interaction.  Thus,  a  weak  shock 
interaction  may  produce  a  limited  separation  bubble  provided  it  does  not 
noticeably  alter  the  global  flow  field.  Shock-boundary  layer  interactions 
are  often  further  subclassified  by  the  character  of  the  shock  and  the 
surface  configuration  such  as:  (1)  incident  oblique  shock  wave,  (2)  com¬ 
pression  corner,  (3)  forward  facing  step,  (4)  rearward  facing  step,  and 
(5)  normal  shock  wave.  Interactions  of  shock  waves  with  laminar  boundary 
layers  are  often  studied  even  though  they  are  rarely  found  in  practical 
applications. 

The  common  occurrence  of  shock-boundary  layer  interactions  and  the 
predominant  nature  of  their  effect  on  flow  phenomena  justify  continued 
effort  in  the  development  of  accurate  prediction  techniques  for  such 
flows.  Green  (Ref  9)  provided  an  excellent  historical  background  (prior 
to  1969)  and  technical  review  of  this  subject  in  a  von  Karnam  Institute 
short  course  on  the  subject  of  shock-boundary  layer  interactions. 
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Schofield  (Ref  23)  also  reviewed  the  subject  and  recommended  an  extensive 
experimental  program  to  advance  existing  understanding  of  such  interactions. 
Analytical,  experimental  and  numerical  approaches  have  been  applied  to  the 
problem  of  developing  prediction  techniques  for  shock-boundary  layer  inter¬ 
actions. 

The  analytical  approach  inherently  requires  that  the  governing  equa¬ 
tions  be  extensively  simplified  before  a  solution  can  be  obtained. 

Analytical  solutions  have  been  developed  for  inviscid  flows,  one-dimen¬ 
sional  viscous  flows,  and  flows  very  close  to  Mach  1.0.  However,  the 
mixture  of  subsonic,  supersonic,  viscous,  and  separation  regions  which 
all  required  simultaneous  treatment  strong  normal-shock  turbulent  boundary 
layer  interaction  results  in  a  fully  elliptic  mathematical  formulation 
which  precludes  analytical  solution  without  extensive  simplification. 

Inger  and  Mason  (Refs  10,  11,  and  15)  have  extensively  studied  shock¬ 
boundary  layer  interactions  analytically.  The  analytical  approach  is  very 
instructive  in  providing  trends,  physical  understanding  of  experimental 
data  and  numerical  results.  Melnick  and  Grossman  (Ref  17)  also  studied 
weak  interactions  with  an  asymptotic/numerical  approach  which  is  a  com¬ 
bination  of  the  analytical  and  numerical  approaches.  However,  it  is  still 
limited  by  simplifying  assumptions  to  the  weak  interaction  case.  At  the 
present  time  no  analytical  method  exists  capable  of  predicting  the  flow 
field  when  separation  occurrs,  even  though  some  analytical  techniques  are 
efficient  for  the  calculation  of  some  critical  characteristics  of  the 
problem  (Ref  18).  These  methods  may  provide  the  conditions  for  Initial 
separation,  but  they  have  little  hope  of  predicting  the  extent  of  separ¬ 
ation  effects  and  would  likely  over-restrict  any  design  process. 
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As  noted  by  Green  (Ref  9)  the  coupled  region  of  the  flow  which  may 
be  regarded  as  "interaction”,  is  more  extensive  when  the  shock  Is  normal. 
The  mixture  of  supersonic  flow,  subsonic  flow,  viscous  boundary  layer, 
vortex  sheet,  strong  normal  shock,  weak  oblique  shocks,  and  possible 
separation  and  reattachment  are  all  exhibited  by  the  normal-shock  inter¬ 
action  problem.  Although  in  the  late  1940's  and  1950's,  extensive 
experimental  effort  (Refs  2,  8,  13  and  19)  concentrated  on  shock  induced 
separation;  the  first  detailed  quantitative  work  for  a  strong  normal-shock 
turbulent  boundary  layer  interaction  was  done  by  Seddon.  Figure  2  illus¬ 
trates  the  strong  normal-shock  turbulent  boundary  layer  interaction  as 
documented  by  Seddon  (Ref  24)  and  has  become  the  accepted  structure  for 
such  interactions.  Recent  experimental  efforts  (Ref  1,  7,  12,  and  16) 
have  added  to  the  base  of  understanding  of  such  interactions.  The  inherent 
problems  of  conducting  experimental  research  at  near  sonic  conditions  have 
limited  the  experimental  results  to  relatively  few  measurements  on  only 
basic  geometries.  The  introduction  of  very  small  probes  into  the  flow  at 
near  sonic  conditions  often  results  in  marked  alteration  of  the  flow. 
Therefore,  the  flow  field  measured  very  frequently  does  not  represent  the 
desired  phenomenon.  Most  experimental  efforts  are  conducted  in  closed 
wind  tunnel  test  sections  where  shock  interactions  with  boundary  layers  on 
side  walls,  floor,  and  ceiling  usually  alter  the  desired  flow  field.  In 
addition  to  the  disturbance  due  to  the  probe,  the  solid  blockage  of  the 
entire  test  apparatus  can  cause  alteration  of  the  flow.  Velocity  measure¬ 
ments  can  be  made  without  introducing  any  apparatus  into  the  test  section 
by  analysis  of  the  Doppler  shift  in  monochromatic  light  when  back-scattered 
from  particles  in  the  flow.  Abbiss  (Ref  1)  and  East  (Ref  7)  have  applied 
such  a  laser  anemometer  to  this  problem  and  their  results  will  serve  as  the 
primary  experimental  basis  for  comparisons  in  the  current  study. 
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Seddon’s  Normal-Shock/Turbulent  Boundary  Layer  Interaction  Model 


The  third  and  currently  most  active  method  of  investigation  is 
numerical  analysis  via  high  speed  computer.  An  extensive  review  of 
finite  difference  methods  in  computational  fluid  dynamics  is  given  by 
Roache  (Ref  21).  Numerical  methods  are  presently  being  applied  to  many 
previously  unsol vable  problems  in  fluid  mechanics.  For  this  effort  a 
finite  difference  scheme  was  chosen  over  finite  element  methods,  because 
the  current  stage  of  finite  element  application  to  aerodynamic  problems 
makes  that  approach  well  beyond  the  scope  of  the  present  study.  . 

Finite  difference  solutions  to  the  Navier-Stokes  equations  by 
asymptotic  convergence  of  the  time  dependent  equations  to  a  steady  state 
value  is  a  well  established  approach  (Ref  3,  5,  6,  14,  16,  20,  21,  and 
25-31).  Finite  difference  methods  are  explicit,  implicit,  or  hybrid 
(a  combination  of  explicit  and  implicit  difference  equations  for  each 
time  step) .  The  implicit  approach  has  an  inherent  numerical  stability 
advantage.  It  allows  larger  time  steps  per  computational  iteration,  but 
usually  requires  a  more  complex  computational  procedure  for  each  time 
step.  An  implicit  approach  similar  to  Roach  (Ref  20)  would  be  very 
attractive  from  a  computational  efficiency  view  point,  but  to  date,  no 
demonstrated  success  of  such  a  method  at  high  Reynolds  numbers  for  a 
transonic  interaction  is  known.  The  hybrid  explicit-implicit  finite 
difference  method  of  MacCormack  (Ref  14) ,  described  in  Appendix  A,  was 
chosen  for  the  current  investigation  because  of  its  proven  success  on 
similar  problems  (Refs  6,  14,  and  25-31)  and  the  fact  that  the  code  was 
readily  available. 
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Investigation  Ob lec fives 

The  primary  objectives  of  this  study  is  three  fold:  (1)  to  further 
validate  the  finite  difference  approach  as  a  strong  tool  applicable  to 
flows  where  sufficient  simplifications  for  analytical  treatment  are  pre¬ 
cluded,  (2)  to  study  the  numerical  simulation  (boundary  location,  boundary 
conditions,  boundary  condition  implementation,  and  mesh  size)  for  possible 
guidance  to  future  applications  of  finite  difference  methods  of  this  type 
to  aerodynamic  problems,  and  (3)  to  study  the  flow  interaction  problem  for 
a  better  understanding  of  the  phenomenon  based  on  the  more  detailed 
information  which  numerical  analysis  results  provide. 

In  summary,  interactions  of  a  normal-shock  with  a  turbulent  boundary 
layer  have  been  calculated  by  the  finite  difference  method  of  MacCormack 
for  various  flow  conditions  and  with  alternate  boundary  condition  implementa¬ 
tion  to  study  both  the  flow  field  and  the  finite  difference  method  and  their 
interrelationship  with  each  other.  Specifically  the  flow  conditions  of 
Abbiss  (Ref  1)  and  East  (Ref  7)  were  duplicated  and  the  computed  solutions 
were  compared  with  experimental  results  for  both  strong  and  weak  interactions. 
Variation  of  boundary  conditions,  implementation  method,  and  mesh  size 
were  performed  to  provide  a  better  understanding  of  the  solution  approach. 
Finally,  an  attempt  was  made  to  understand  some  of  the  three-dimensional 
effects  in  the  experimental  data  used  for  comparison. 
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II 


This  numerical  investigation  of  normal-shock  turbulent  boundary 
layer  interactions  was  accomplished  using  the  same  basic  code  as  Shea 
(Ref  28) .  The  scheme  is  basically  the  hybrid  implicit-explicit  finite 
difference  method  of  MacCormack  (Ref  14) ,  applied  to  the  time-dependent 
two-dimensional  Navier-Stokes  equations  (see  Appendix  A) .  Boundary 
conditions  and  their  implementation  constituted  a  major  portion  of  the 
effort  in  this  investigation.  First  a  base  line  case  was  established, 
then  the  investigation  proceeded  to  include  varying  Mach  number,  boun¬ 
dary  placement,  and  mesh  size  refinement  with  a  common  boundary  condition 
formulation.  Finally,  three-dimensional  effects  were  addressed  to  a 
limited  extent  in  an  attempt  to  improve  the  comparisons  of  computation 
with  experiment. 

Base  Line  Case 

The  first  task  of  the  investigation  was  to  obtain  a  converged  solu¬ 
tion  for  the  reference  or  base  line  case.  The  flow  conditions  of  Ref  1, 
freestream  Mach  number  of  1.51  and  Reynolds  number  of  30  x  10^  per  ft,  were 
chosen.  The  approach  was  to  establish  a  single  reference  case  and  make 
primary  comparisons  with  experiment  for  that  base  line  case.  Then  related 
studies,  deviating  from  the  base  line  case,  were  compared  with  the  base 
line  results  to  evaluate  the  effect  of  the  change.  A  review  of  the 
criterion  for  the  Ref  1  data  being  chosen  as  the  base  line  case  was  pre¬ 
sented.  Then  a  review  of  the  computational  approach  used  to  compute  the 
base  line  case  was  followed  by  an  outline  of  the  other  studies  and  com¬ 
parisons  to  the  base  line  case. 

(  The  data  of  Seddon  (Ref  24),  which  were  obtained  on  a  flat  plate 

mounted  in  a  transonic  test  section,  was  the  basis  for  the  experimental 


comparisons  of  Shea.  The  shock  was  generated  by  a  second  plate  with  a 
flap  mounted  above  and  downstream  of  the  flat  plate  leading  edge  (see 
Figure  3a).  The  effect  of  the  shock  generator  on  the  downstream  portion 
of  the  interaction  region  and  the  unknown  boundary  layer  characteristics 
on  the  side  walls  probably  affect  the  flow  field  in  ways  that  are  difficult 
to  predict.  The  flat  plate  was  relatively  short  with  the  primary  measuring 
station  approximately  9  inches  from  the  leading  edge.  Thus,  the  boun¬ 
dary  layer  thickness  at  the  interaction  was  approximately  0.16  inches 
(Ref  28:5-6).  On  the  other  hand,  the  Abbiss  data  (Ref  1)  was  collected  by 
a  configuration  where  the  shock  was  generated  across  the  entire  test  sec¬ 
tion  by  a  variable  nozzle  in  the  tunnel  diffuser  (see  Figure  3b).  The 
interaction  generated  by  this  configuration  has  the  following  inherent 
advantages  over  the  Seddon  data: 

a)  The  boundary  layer  is  approximately  1.6  inches  thick,  10 
times  the  Seddon  test,  at  the  interaction  zone. 

b)  The  shock  generator  is  far  downstream  and  should  not  effect 
the  Interaction  region. 

c)  The  boundary  layer  should  be  very  similar  on  all  four  walls, 
therefore,  the  effects  should  be  symmetrical  and  more  pre¬ 
dictable. 

Data  produced  by  tests  (Ref  4)  incorporating  a  wall  bump  to  produce 
the  shock  wave  (Figure  3c)  were  not  selected.  The  unsymmetric  test  sec¬ 
tion  configuration  and  pressure  gradients  downstream  of  the  shock  would 
introduce  undesired  complications  to  the  flow  field. 

The  laser  anemometer,  which  was  the  primary  means  of  data  acquisi¬ 
tion  for  both  Abbiss  and  East,  does  not  introduce  any  probes  or  support 
mechanicisms  into  the  test  section.  Therefore,  the  likelihood  of  alter¬ 
ing  the  flow  field  by  the  measuring  devices  was  almost  totally  eliminated. 

The  Mach  1.51  case  was  chosen  as  the  base  line  case  because  no  pre¬ 
vious  comparisons  to  these  data  are  known,  the  test  approach  provides  a 
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flow  field  which  should  be  more  symmetric  than  the  other  known  data, 
there  are  more  data  available  for  this  case  than  for  any  others  with  this 
test  configuration,  and  the  data  acquisition  approach  is  expected  to 
have  no  effect  on  the  flow  field.  The  Mach  1.51  flow  is  very  comparable 
to  the  results  of  Seddon  and  Shea  (M  =  1.47)  which  can  be  used  to  provide 
additional  information  and  understanding.  With  the  flow  conditions  and 
finite  difference  approach  established,  attention  was  focused  on  establish¬ 
ing  the  computational  mesh  structure  and  boundary  conditions  and. then 
obtaining  a  converged  solution  for  the  base  line  conditions. 

The  computational  domain  was  extended  in  the  y  direction  to  1.5  ft 
from  the  wall  compared  to  1.33  ft  by  Shea  (Ref  29).  This  extended  outer 
edge  of  the  computational  domain  is  geometrically  at  the  same  distance 
from  the  wall  as  the  wind  tunnel  center  line.  Thus,  the  computational 
domain  is  bounded  by  a  plane  of  symmetry  and  a  flat  wall.  The  boundary 
conditions  for  the  base  line  case  are  presented  here  for  the  five  regions 
shown  in  Figure  4. 

Wall  Boundary.  The  wall  boundary  is  prescribed  after  Shea  (Ref  28) . 
The  boundary  conditions  are: 

No  slip  wall 

u  -  v  =  0  11-1 

Adiabiatic  wall 

(T)y  =  o  H-2 

The  compatible  condition  of  the  normal  (y)  momentum  equation  at  the  wall 
yields  (Ref  26)  a  balance  of  the  pressure  and  shear  gradients  at  the  sur¬ 
face. 

<p)y  »  <b(u)  )  +  (4/3  p(v)  )  -  (2/3  p(u)  )  H-3 

7  y  x  '  y  '  X  y 
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Subsonic  Center  Line  Boundary 
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Boundary  Condition  Region  Definitions 


The  boundary  conditions  were  implemented  by  a  reflection  scheme  where 
the  first  two  rows  (J  =  1  and  J  s  2)  of  mesh  points  are  equal  equidis¬ 
tant  from  the  wall  on  opposite  sides.  The  J  «  1  row  is,  therefore, 
positioned  outside  the  flow  field  and  given  values  which  would  result  in 
the  desired  values  at  the  surface  by  linear  interpolation.  The  resulting 
finite  difference  expressions  are  presented  in  Appendix  B. 

Inflow  Boundary.  The  inflow  boundary  values  are  from  a  finite  dif¬ 
ference  boundary  layer  code  developed  by  Shang  and  documented  by .Beauregard 
(Ref  5).  The  values  for  u,  v,  p,  and  p  for  all  1=1  (inflow  boundary) 
mesh  points  were  interpolated  from  the  output  of  the  boundary  layer  code. 
These  values  are  held  fixed. 

Supersonic  Center  Line  Boundary.  The  supersonic  portion  of  the 
center  line  boundary  is  also  in  a  fully  supersonic  region  and  the  input 
is  held  constant  with  the  same  value  as  the  center  line  point  on  the 
inflow  boundary.  No  changes  were  made  relative  to  Shea  for  the  wall, 
inflow,  or  supersonic  center  line  boundary  treatment. 

Subsonic  Center  Line  Boundary.  The  upper  subsonic  boundary  based  on 
the  new  center  line  condition  results  in  a  formulation  quite  different 
from  that  of  Shea.  The  symmetry  condition  at  the  upper  boundary  provides 
a  number  of  possible  boundary  conditions.  All  derivatives  of  state  pro¬ 
perties  with  respect  to  spatial  variable  perpendicular  to  the  line  of 
symmetry  (y)  are  zero. 

(P)y  =  <P)y  =  <ei>y  -  <et>y  =  (T)y  =  0  II-4 

The  variation  in  y  of  all  vector  quantities  are  odd  functions  for  vector 
components  perpendicular  to  the  center  line: 

v(y)  -  -v(-y)  11-5 

where  y  is  measured  from  the  center  line.  Similarly,  the  variation  in 


y  of  all  vector  quantities  are  even  functions  for  vector  components 
parallel  to  the  center  line; 

u(y)  *  u(-y) 


thus 


(u)  v  =  0 

y 


11-6 

II-7 


If  the  flow  at  the  center  line  is  also  assumed  to  be  lnvlscid  and 
adiabatic,  the  steady  state  solution  can  be  shown  (see  Appendix  B)  to 
require  that 


(et)x  -  0  '  II-8 

Several  boundary  condition  formulations  are  therefore  possible  for  the 
center  line  boundary.  The  base  line  case  was  computed  with 

v  =  (u)y  =  (p)y  =  (p)y  =  0  11-9 

The  implementation  of  the  center  line  boundary  derivative  quantities 
is  by  a  two-point  first-order  extrapolation.  Letting  Z  be  a  typical  var¬ 
iable  (see  Appendix  B) 

Z(I,  JL)  =  Z(I.  JLMI)  +  Z(I.  JLM2)  11-10 

2 

Which  gives  a  value  of  u,  p,  and  p  on  the  subsonic  center  line  which  is 
an  average  of  the  two  points  just  below  it.  Alternate  formulations  of 
the  center  line  boundary  conditions  and  their  implementation  will  be  dis¬ 
cussed  in  Section  III. 

Outflow  Boundary.  The  outflow  boundary  has  been  found  to  be  rather 
important.  The  base  line  case  was  computed  using  an  assumption  that  the 
outflow  boundary  is  far  enough  downstream  that  the  primitive  variables 
(p,  u,  v,  and  p)  are  nearly  invariant.  The  no-change  condition  is  commonly 
used  for  MacCormack's  scheme  to  minimize  the  upstream  influence  (Ref  14,  26). 


(p)x  ■  (u)x  -  (v)x  *  (p)x  •»  0  11-11 

These  conditions  were  implemented  by  a  one-point  first-order  extrapolation 
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Again,  letting  Z  be  a  typical  variable 

Z(1L,  J)  =  Z(ILM1 ,  J)  II"12 

for  all  J.  Discussion  of  alternate  formulations  of  the  outflow  boundary 
is  provided  in  Section  111. 

Other  Studies 

Two  other  investigations  were  conducted.  They  were  variation  of 
outflow  boundary  position  and  variation  of  free  stream  Mach  number. 

Outflow  Boundary  Position.  After  the  base  line  case  was  established, 
variation  of  the  outflow  boundary  position  was  investigated.  The  outflow 
boundary  was  located  at  three  alternate  positions.  The  base  line  location 
of  the  outflow  boundary  is  approximately  18  times  the  undisturbed  boundary 
layer  thickness  downstream  of  the  shock  location.  The  first  alternate 
position  was  half  way  between  the  base  line  position  and  the  shock  wave 
or  approximately  9  boundary  layer  thicknesses  downstream  from  the  shock. 
Next,  the  outflow  boundary  was  placed  midway  between  the  first  two  cases 
or  approximately  13  to  14  boundary  layer  thicknesses  behind  the  shock 
location.  The  streamwise  mesh  spacing  was  the  same  for  both  of  these  for¬ 
ward  locations  as  it  was  for  the  base  line  case.  Because  of  convergence 
problems  with  both  of  the  cases  for  outflow  boundaries  forward  of  the  base 
line  case,  a  third  case  was  computed. 

The  third  alternate  outflow  boundary  position  was  approximately  29 
boundary  layer  thicknesses  from  the  shock  location.  This  location  was 
obtained  by  increasing  the  streamwise  mesh  spacing  by  50  percent.  No 
attempt  was  made  to  run  a  smaller  mesh  size  because  of  time  constraints 
and  computer  memory  limitations  with  the  larger  number  of  mesh  points 
required . 
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Mach  Number  Variation 

In  addition  to  the  base  line  case  (M  **  1.51)  cases  lor  >1  *=1.3  and 

00  00 

=  1.4  forward  of  the  normal  shock  were  also  computed.  These  cases 
were  calculated  using  the  same  boundary  conditions  as  for  the  base  line 
case.  The  results  of  these  simulations  are  compared  with  the  experimental 
results  of  East  (Ref  7)  in  the  next  section. 


Ill  Discussion  and  Results 

The  computed  base  line  case,  «  1.51,  was  compared  with  the 
experimental  data  of  Abbiss  (Ref  1).  Several  alternate  boundary  condition 
formulations  were  computed  and  compared  with  both  base  line  computation 
and  the  Abbiss  experimental  data.  Computations  were  also  performed  for 
=  1.30  and  =  1.40  with  the  same  boundary  condition  formulation  as 
the  base  line  case,  and  these  results  were  compared  with  the  experimental 
results  of  East  (Ref  7) .  The  East  results  were  obtained  in  a  subsequent 
test  in  the  same  facility  as  that  of  Abbiss  with  only  minor  modifications 
to  the  experimental  approach.  Since  separation  is  inherently  a  three- 
dimensional  phenomenon,  three-dimensional  effects  were  discussed.  Results 
from  an  attempt  to  correct  the  two-dimensional  numerical  method  for  the 
most  significant  three-dimensional  effects  were  presented  and  compared 
with  the  computations  of  Shea  (Ref  29). 

Identification  of  Results.  An  indentif ication  code  was  assigned  to 
each  case  since  several  different  cases  were  computed.  Table  1  gives  a 
complete  list  of  indentif ication  codes  and  a  corresponding  description 
of  each  series  of  computations.  The  alphabetic  part  of  the  code  identifies 
the  flow  conditions  and  boundary  condition  formulation  with  respect  to  the 
base  line  case.  The  numerical  portion  of  the  identification  indicates 
subsequent  iteration  of  the  solution.  The  higher  numbers  are  later  runs. 
Typically,  the  solutions  were  advanced  120  time  steps  for  a  total  elapsed 
time  of  approximately  0.002  to  0.0025  seconds  between  data  sets.  Table  1 
is  provided  as  a  guide  to  the  identification  of  the  data  source  for  data 
plots  in  this  section. 

Solution  Convergence.  A  convergence  criterion  was  established  and 
applied  as  follows.  The  wall  surface  pressure  distribution  was  plotted 
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Table  I 


Computation  Code  ID  Key 


Description  and  Remarks 

Base 

Line  Case,  See  Section  II  For  Description 

CLBC1 

< 

II 

•O 

li 

/"N 

P 

II 

li 

o 

Type 

1 

OFBC2 

(P)x  =  (Pa)x  =  (Pv)x  =  (et)x  =  0 

Type 

2 

CLBC 

V  =  <P)y  =  (u)y  =  (p)y  =  0 

Type 

3 

OFBC 

(u)  =  (v)  =  (p)  =  (p)  =  0 

Type 

4 

XX  XX  XX  XX 

CLBC 

All  Values  Fixed  at  Normal-Shock 

Values 

OFBC 

Same  as  H-Series 

CLBC 

Same  as  H-Series 

OFBC 

Same  as  O-Sories 

CLBC 

V  =  (p)  =  (pu)  =  (H  )  =  0 

K  y  y  t  x 

ofbc 

Same  as  0-Series 

CLBC 

Same  as  H  Except  Type  5 

OFBC 

Same  as  0-Series 

CLBC  V  =  (p)y  =  (Ht)x  -  (pu)x  =  0 


OFBC  Same  as  0-Series 


All  Boundary  Conditions  Same  as  0-Series 
IL  =  40  (IL  -  64  for  all  previous  cases) 


Table  I  (cont'd) 


R 

All  Boundary  Conditions  Same  as  0-Series 

IL  =  52 

S 

Same  as  0-Series  Except  Dx  =  0.778  (50%  increase  from 
value  for  all  other  Series) 

Q 

All  Boundary  Condition  Same  as  0-Series 

Mro  =  1.30 

V 

Boundary  Conditions  Same  as  0-Series 

M  =  1.40 

OO 

B 

See  Reference  28  (B-Series  Constitutes  Ref  29) 

I 

All  Boundary  Conditions  Same  as  0-Series 

JL  =  22  (JL  =  32  All  other  runs  except  U-Series) 

T 

Same  as  0-Series 

Except: 

JL 

Y  Pu(J)  DY(J) 

J=1  1  =  1 

JL 

=  y  pu(J)  DY(J) 

J.1  1 2 3  "  1L 

U 

Same  as  I-Series 

1.  Subsonic  Center  Line  Boundary  Conditions 

2.  Outflow  Boundary  Conditions 

3.  See  Table  2  for  description  of  difference  equation  types 


Table  II 


every  30  time  steps.  When  the  wall  pressure  ratios  remained  constant, 
did  not  change  by  more  than  0.01  (0.3  -*•  0.8%)  in  a  characteristic  time 
(approximately  180  time  steps),  convergence  was  said  to  be  established. 

With  an  expected  accuracy  of  the  experimental  data  no  better  than  1%  and 
the  expected  comparison  between  computation  and  experiment  approximately 
5  to  10%,  the  0.01  criterion  represents  an  attainable  resolution  which 
should  not  introduce  scatter  into  the  results  of  this  study.  Early  in 
the  investigation,  complete  longitudinal  distributions  were  plotted  for 
sequential  data  sets.  It  was  concluded  that  the  pressure  variation  of 
selected  locations  plotted  verses  time  gave  a  good  representation  of  the 
steady  state  asymptote.  Figure  5  is  the  convergence  criterion  plot  for 
the  base  line  (0  Series)  results.  It  is  a  typical  result  and  any  devia¬ 
tion  from  this  type  of  plot  will  be  discussed. 

The  results  presented  herein  were  obtained  using  the  (Wright-Patterson 
AFB,  ASD)  Cyber  175  computer.  Computations  required  approximately  15 
minutes  per  characteristic  time.  Thus,  60  to  75  minutes  were  required  to 
obtain  a  converged  solution  with  a  boundary  condition  formulation  similar 
to  the  base  line  case. 

Base  Line  Case  Results 

The  numerical  and  experimental  (Ref  1)  results  show  qualitative 
agreement  as  can  be  seen  in  Figure  6.  The  lambda  shock  structure  can 
clearly  be  seen  in  the  computed  Mach  contours  of  Figure  6b.  The  boundary 
layer  exhibits  rapid  growth  and  a  series  of  compression  waves  are  gener¬ 
ated  in  the  supersonic  flow.  Although  these  compression  waves  do  not 
form  a  distinct  bifurcation  point,  they  are  collimated  into  a  single 
longitudinal  mesh  space  at  approximately  the  same  height  as  the  experi¬ 
mentally  determined  bifurcation  point.  The  shear  layer  which  is  shed 
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from  the  bifurcation  point  is  not  observed  in  the  numerical  results 
because  the  mesh  spacing  is  too  coarse  in  that  region. 

The  Abbiss  results  do  not  provide  the  details  of  the  boundary  layer 
in  the  interaction  region  because  of  measurement  difficulties  close  to 
the  wind  tunnel  wall.  Seddon  (Ref  24),  however,  shows  that  for  =  1.47 
the  separation  bubble  length  is  approximately  12  times  the  undisturbed 
boundary  layer  thickness  which  agrees  with  the  computed  results.  The 
computations  and  Seddon  both  show  a  smooth  pressure  distribution  in  the 
reattachment  region.  The  maximum  separation  bubble  height  and  the  loca¬ 
tion  of  the  maximum  height  of  separation  in  the  Seddon  data  and  the  present 
numerical  results  are  quite  different.  Seddon' s  results  show  a  maximum 
separation  bubble  height  of  about  half  the  initial  boundary  layer  thickness 
compared  to  approximately  0.1  the  initial  thickness  height  for  the  com¬ 
putations.  Seddon  measured  the  maximum  height  at  70%  of  the  bubble  length 
compared  to  a  40%  location  for  the  computed  results.  The  poor  agreement 
should  not  be  unexpected.  Since  the  separation  bubble  size  and  structure 
is  likely  a  strong  function  of  the  post  shock  pressure  distribution,  and 
the  post  shock  pressure  distribution  is  probably  affected  by  the  differences 
between  the  Seddon  experimental  arrangement  and  the  numerically  simulated 
problem  which  is  based  on  the  Abbiss  test  configuration.  In  addition,  the 
turbulence  model  was  an  algebraic  mixing  length  model  developed  for  primary 
application  to  attached  turbulent  boundary  layers  in  equilibrium.  Although 
the  use  of  this  model  (Ref  28)  and  similar  models  (Refs  5,  6,  14,  16,  25-31) 
is  common  practice,  they  should  not  be  expected  to  model  the  detailed 
separation  structure.  Viegas  et  al  (Refs  6,  30,  31)  have  studied  more 
complex  turbulence  models  extensively  with  mixed  results.  Further,  separ¬ 
ation  and  turbulence  are  inherently  a  three-dimensional  phenomena  and 
detailed  modeling  cannot  be  accomplished  by  a  two-dimensional  analysis. 
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The  most  obvious  disagreement  which  can  be  observed  in  figure  6 
is  the  total  lack  of  a  post  shock  supersonic  region,  commonly  called  the 
supersonic  tongue.  Three  possible  explanations  for  this  discrepancy 
are  offered.  First,  the  mesh  structure  is  quite  coarse  in  the  region 
where  the  supersonic  tongue  is  thought  to  exist  and  the  resolution  may 
not  be  sufficient  to  exhibit  the  required  detail  to  describe  the  super¬ 
sonic  region.  Second,  the  supersonic  region  is  quite  possibly  a  product 
of  a  three-dimensional  blockage  effect  produced  by  the  rapid  boundary 
layer  growth  on  all  four  walls  of  the  post  shock  region  of  the  interact¬ 
ion.  Finally,  the  data  reduction  technique  assumed  constant  static  pres¬ 
sure  in  the  boundary  layer  from  which  local  Mach  numbers  were  calculated. 

The  computed  local  flow  variables  in  the  post  shock  boundary  layer  region 
are  not  constant.  The  differences  are  of  sufficient  magnitude  to  produce 
local  Mach  numbers  of  slightly  less  than  M  «*  1.0  in  the  post  shock  outer 
boundary  layer  edge  region.  It  should  also  be  noted  that  measurements  in 
this  region  required  that  the  optics  pass  through  the  strong  interaction 
region  of  the  tunnel  side  walls  and  that  the  Seddon  model  was  used  by  Abbiss 
as  a  guide  for  his  experimental  effort.  Also,  no  evidence  of  post  shock 
supersonic  flow  can  be  discerned  from  the  Schlieren  data  in  the  Seddon 
results. 

A  more  quantitative  comparison  of  the  base  line  computation  with 
the  Abbiss  data  is  presented  in  Figures  7  and  8.  As  will  be  discussed 
further  in  a  later  section,  the  shock  could  not  be  positioned  exactly  in 
a  given  longitudinal  location  by  the  shock  Implementation  procedure  used 
for  the  computations.  Although  the  relative  location  of  the  shock  varied 
less  than  one  undisturbed  boundary  layer  thickness  (not  enough  to  affect 
the  interaction)  from  the  established  reference  location  for  most  cases 
computed,  it  would  cause  unrealistc  scatter  in  the  results  if  not  resolved 
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Figure  7.  Base  Line  Case  Experimental  and  Numerical  Longitudinal 
Total  Velocity  Profiles 
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Numerical  Transverse  Total  Velocity  Profiles 


prior  to  plotting.  For  one  profile  where  the  shock  jump  was  distinctive 
with  respect  to  x-location  (y/$1  *  4.954)  the  effective  x-location  was 
adjusted  for  all  computed  and  measured  results.  This  relative  location 
was  then  held  fixed  for  all  subsequent  plotting. 

The  comparison  of  longitudinal  profiles  (Fig  7)  shows  practically 
all  points  are  within  10%  of  the  measured  values.  The  agreement  is 
excellent  for  the  y/S^  =  0.944  curve  which  is  closest  to  the  wall.  The 
superior  agreement  for  this  profile  is  attributed  to  the  refined  .mesh  in 
this  region  and  the  reduction  of  three-dimensional  effects  due  to  the 
proximity  of  the  wall. 

The  experimental  velocities  are  consistently  higher  in  the  subsonic 
regions  of  the  flow.  The  velocity  decrease  across  the  shock  is  approxi¬ 
mately  100  ft/sec  less  for  the  experimental  results  relative  to  the 
calculated  values.  The  shock  jump  when  predicted  by  the  one-dimensional 
Rankine-Hugoniot  relations  yields  a  post  shock  value  of  766  ft/sec.  If 
instead  of  a  normal-shock  an  81  deg  oblique  shock  jump  is  computed,  the 
results  agree  well  with  the  measured  data  and  an  85  deg  oblique  shock  jump 
agrees  with  the  computed  results.  Therefore,  the  differences  noted 
between  the  measured  and  computed  results  can  be  accounted  for  by  small 
differences  in  the  ..hock  angles  relative  to  the  local  flow  in  the  super¬ 
sonic  regions.  This  explanation  would  also  account  for  the  much  better 
agreement  in  the  single  subsonic  profile  (y/ 6^  =  0.944). 

The  subsonic  region  of  the  profiles  shows  excellent  agreement  in 
shape.  Relative  to  the  calculated  values,  higher  experimental  velocities 
are  found  in  the  supersonic  regions  and  lower  experimental  velocities  are 
in  the  subsonic  region.  The  consistently  higher  supersonic  velocities  and 
lower  subsonic  velocities  suggest  the  effect  of  an  area  restriction  on 
both  regions  of  the  flow.  It  is  concluded  that  the  predominant  difference 
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between  the  experimental  results  of  Abbiss  and  the  computed  results  is 
the  three-dimensional  area  restriction  imposed  or  the  flow  by  the  rapid 
expansion  of  the  boundary  layer  displacement  thickness  on  the  side  walls 
of  the  tunnel.  This  three-dimensional  effect  is  not  modeled  by  the  two- 
dimensional  computations  reported  herein. 

The  vertical  profiles  of  Figure  8  present  the  same  data  in  a  dif¬ 
ferent  format  and  support  the  same  conclusions.  Once  again  the  vertical 
profiles  show  qualitative  agreement  in  shape.  The  additional  blockage 
due  to  three-dimensional  effects  is  also  a  likely  cause  of  the  higher 
shock  angles  relative  to  the  flow  direction  in  the  experimentally  measured 
results.  Shea  (Ref  28)  argued  that  the  Seddon  results  exhibited  strong 
three-dimensional  effects.  This  claim  was  supported  by  Winter  (Ref  32), 
who  indicated  that  the  interaction  region  for  both  the  Abbiss  and  East 
data  are  strongly  three-dimensional..  Winter's  statement  was  based 
on  oil  flow  and  pressure  data  obtained  in  conjunction  with  the  experi¬ 
mental  efforts  reported  by  Abbiss  and  East.  While  the  computed  results 
are  different  from  the  experimental  data,  all  of  the  differences  are  con¬ 
sistent  with  the  expected  difference  between  the  two-dimensional  formu¬ 
lation  of  the  base  line  case  and  the  measured  data  which  exhibit  known 
three-dimensional  effects.  An  attempt  to  account  for  these  three-dimen¬ 
sional  effects  will  be  discussed  later. 

Boundary  Conditions 

Outflow  Boundary  Conditions.  The  early  phase  of  the  investigation  con¬ 
centrated  almost  completely  on  the  evaluation  of  alternate  subsonic 
center  line  boundary  treatment  and  the  outflow  boundary  condition  was 

(V)X  X  "  (u)x  X  -  <p>x  X  "  (f°x  X  "  0  m“l 
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Figure  9a  shows  a  Mach  contour  plot  of  the  final  H-series  computation. 

The  particular  center  line  boundary  treatment  for  this  series  is  defined 
in  Table  I.  for  completeness  but  is  not  pertinent  here.  The  plot  is  pre¬ 
sented  as  typical  of  all  results  with  the  second-order  outflow  boundary 
condition.  The  outflow  pressure  continued  to  rise  until  the  shock  wave 
anA  the  interaction  was  crowded  forward  into  the  fixed  inflow  boundary 
and  the  outflow  decreased  to  M  =  0.25  which  was  also  unreasonable.  The 
outflow  conditions  of  Eqn  III-l  were  recommended  in  Roache  (Ref  21) ,  but 
have  been  determined  to  be  unacceptable  for  a  subsonic  outflow  boundary, 
at  least  in  conjunction  with  the  MacCormack  code  used  herein.  In  an  attempt 
to  understand  the  problem,  the  subsonic  center  line  condition  was  fixed 
with  the  Rankine-Hugoniot  normal  shock  values  at  all  post  shock  mesh  points. 
The  J-series  (Fig  9b)  Mach  contours  resulted.  A  comparison,  of  the  H  and  J 
series  computations,  indicated  that  the  result  of  changing  the  center  line 
boundary  treatment  was  to  stablize  the  outflow  pressure  and  velocity. 
Subsequently,  the  outflow  boundary  condition  was  changed  to  the  first 
derivative  condition. 

(v)x  =  (P)x  =  (u)x  =(p)x  =  0  III-2 

No  claim  is  intended  that  the  boundary  conditions  of  Eqn  III-2  are  the 
product  of  an  exhaustive  study  or  that  they  are  the  best  treatment  pos¬ 
sible.  These  conditions  did  provide  a  series  of  solutions  to  converge 
numerically  to  values  which  agree  quite  well  with  experimental  results. 

As  was  described  in  Section  II,  the  outflow  boundary  conditions  of  Eqn 
HI-2  were  implemented  by  assigning  the  values  on  the  outflow  boundary 
the  same  values  as  those  of  the  adjacent  points  just  upstream.  All  com¬ 
putations  presented  herein  except  the  H,  J,  T,  and  U  series  were  computed 
with  this  simple  first-order  scheme  applied  to  the  primitive  variables 
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u,  v,  p  and  p.  Series  T  and  U  formulation  and  results  will  be  discussed 
in  the  three-dimensional  effects  section. 

Center  Line  Boundary  Conditions.  Table  1,  page  19,  outlines  in  detail 
the  center  line  boundary  conditions  of  Figure  10.  The  subsonic  center 
line  boundary  was  treated  by:  a  first-order  single-point  extrapolation 
(K-series) ,  a  two-point  second-order  scheme  (M-series) ,  a  constant  total 
enthalpy  and  momentum  approach  (N-series)  and  the  two-point  first-order 
averaging  scheme  of  the  base  line  (0-series).  The  results  (Fig  IT)  show 
no  discernible  difference  with  the  possible  exception  of  the  base  line  case 
for  the  x/6^  =  79.5  profile.  The  conclusion  is  that  any  compatible  formula¬ 
tion  of  the  center  line  boundary  will  produce  a  reasonable  result  if  com¬ 
bined  with  an  appropriate  outflow  boundary  treatment. 

The  base  line  results  for  the  x/6^  =  79.5  profile  exhibit  the  only 
measurable  difference  from  the  other  computations  in  Figure  10.  This 
deviation  of  a  single  profile  just  aft  of  the  shock  location  can  be  traced 
to  the  shock  implementation  procedure.  The  normal  shock  jump  at  the  center 
line  is  forced  to  meet  the  Rankine-Hugoniot  relations  between  the  supersonic 
and  subsonic  regions  of  the  center  line  boundary.  Figure  11  presents  two 
sequential  Mach  contour  plots  to  examine  the  shock  jump  implementation 
at  the  center  line  boundary.  The  K-series  was  selected  as  typical  of  the 
early  shock  implementation  scheme.  The  K7  contours  in  Figure  11a  show 
only  the  M  =  0.8  contour  embedding  the  single  point,  fixed  by  the  Rankine- 
Hugoniot  relations,  just  downstream  of  the  shock  on  the  center  line  boun¬ 
dary.  Less  than  1.5  characteristic  times  later,  the  K9  contours  of  Figure 
lib  show  both  the  M  ■  0.8  and  0.9  contours  have  passed  the  fixed  point  and  a 
measurable  downstream  movement  in  the  entire  shock  structure  can  be  observed. 
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The  wall  pressure  ratios  on  the  convergence  plot  also  revealed  a  slow 
steady  decrease  as  the  shock  structure  moves.  If  allowed  to  continue, 
the  entire  shock  will  slowly  bleed  past  the  fixed  point.  Then  the  uncon¬ 
strained  shock  wave  will  almost  immediately  pass  completely  out  of  the 
computational  domain.  In  order  to  prevent  this  and  to  hold  fixed  the 
shock  jump  on  the  upper  boundary,  five  consecutive  points  downstream  of 
the  desired  shock  location  were  assigned  values  computed  by  the  Rankine- 
Hugoniot  relations. 

Several  other  approaches  were  tried.  The  shock  was  allowed  the 
freedom  to  spread  over  3  or  4  mesh  points  but  it  consistantly  formed 
either  at  the  forward  or  rear  mesh  space  of  the  allowable  range.  The 
procedure  to  fix  five  center  line  boundary  points  with  the  Rankine- 
Hugoniot  post  shock  values  was  established  for  the  base  line  case  and 
used  in  all  subsequent  computations,  The  x/6^  =  77.2  center  line  point 
for  the  base  line  case  is  in  this  five-point  range  and  therefore,  is 
artifically  held  at  a  lower  value  of  local  velocity  than  would  otherwise 
be  predicted.  This  procedure  did  allow  the  consistent  convergence  of 
the  numerical  method  and  only  local  effect  are  exhibited  by  Figure  10. 

Outflow  Boundary  Location.  The  first-order  boundary  condition  of 
Eqn  III-2  was  used  with  a  common  set  of  center  line  boundary  conditions 
(same  as  base  line  case)  to  compute  results  for  three  alternate  outflow 
boundary  locations.  The  most  forward  location  for  the  outflow  boundary 
was  approximately  9  times  the  undisturbed  boundary  layer  thickness  down¬ 
stream  of  the  shock  wave.  This  location  is  denoted  as  the  P-series  and 
was  computed  at  only  40  (IL=40)  streamwise  mesh  points.  The  second 
position  (denoted  the  R-series)  placed  the  outflow  boundary  midway 
between  the  most  forward  position  and  the  base  line  location  which  was 
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approximately  18  times  the  initial  boundary  layer  thickness  from  the  shock 
wave.  Finally,  the  streamwise  mesh  spacing  was  increased  by  50%,  which 
allowed  placement  of  the  outflow  boundary  at  approximately  29  times  the 
undisturbed  boundary  thickness  downstream  of  the  shock  wave  without 
increasing  the  number  of  mesh  points.  Figure  12  provides  a  composite 
plot  of  the  velocity  profiles  for  all  of  these  outflow  boundary  locations. 
Only  the  second-order  outflow  boundary  treatment  shows  measurable  devia¬ 
tion  from  the  base  line  results.  Since  neither  changes  in  outflow  boundary 
location  nor  increased  streamwise  mesh  spacing  produce  a  measurable  effect 
on  the  solution,  it  is  concluded  that  the  streamwise  mesh  is  small  enough 
and  the  outflow  boundary  is  sufficiently  far  downstream  that  neither  is  a 
strong  influence  on  the  solution  in  the  interaction  region.  A  closer 
examination  of  the  reattachment  point  and  the  wall  shear  stress  values 
revealed  that  both  were  measurably  altered  by  the  two  forward  locations. 
Therefore,  the  outflow  boundary  should  not  be  moved  forward  of  the  base 
line  case  without  examination  of  potential  influence  on  the  computed  results. 

Mach  Number  Variation 

The  same  boundary  condition  formulation  as  the  base  line  case  was 
used  to  compute  the  =  1.30  and  =  1.40  cases.  These  cases  were 
computed  in  a  production  manner  and  are  compared  to  the  East  (Ref  7)  results 
in  Figure  13.  The  comparisons  show  better  agreement  than  the  Mach  1.51  case. 
The  vertical  profiles  of  Figure  13a  show  agreement  within  5%  to  7%  and 
qualitative  agreement  in  trends  and  shape  for  the  corresponding  curves. 

Once  again  the  higher  subsonic  values  and  smaller  shock  jumps  in  the  experi¬ 
mental  data  indicate  measurable  three-dimensional  effects.  The  results  are 
consistent  with  the  separation  criterion  of  Schlichting  (Ref  22:365)  which 
states  that  M^  1.30  is  required  to  promote  separation.  The  Mw  =  1.40 
case  shows  a  small  separation  bubble  of  length  approximately  4  times  the 
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undisturbed  boundary  layer  thickness.  The  Mre  »  1.30  computations  had  an 
even  shorter  bubble  which  was  less  than  3  times  the  undisturbed  boundary 
layer  thickness  length  with  a  computed  maximum  bubble  height  of  0.002  ft 
(less  than  0.02  6^)  which  would  be  extremely  difficult  to  observe  and 
was  obviously  well  into  the  weak  interaction  classification. 

Three-Dimensional  Effects 

As  was  previously  stated,  separation  is  not  two-dimensional  and 
therefore  inherently  exhibits  three-dimensional  effects.  The  transonic 
character  of  the  flow  field  studied  caused  disturbances  to  be  perturbed 
across  the  entire  duct  by  relatively  small  displacements  on  any  of  the 
boundaries.  Shea  attempted  to  account  for  the  predominant  three-dimen¬ 
sional  effect  by  a  boundary  layer  blockage  computation  which  was  used  to 
determine  the  longitudinal  pu  distribution  on  the  subsonic  region  of  the 
upper  boundary.  This  procedure  was  outlined  in  Reference  28  and  was 
incorporated  xn  the  Shea  computations  of  the  base  line  flow  conditions 
(Ref  29).  The  Shea  results  are  denoted  as  the  B-series  on  Figure  14.  A 
close  examination  of  Figure  14  revealed  that  no  discernible  improvement 
over  the  present  base  line  case  could  be  seen  from  the  Shea  upper  boundary 
treatment. 

The  idea  of  accounting  for  total  boundary  layer  blockage  was  approached 
differently  in  the  current  effort.  A  major  factor  in  any  duct  flow  is  to 
model  the  correct  area  distribution.  If  the  area  distribution  is  duplicated 
for  different  geometric  shapes,  very  similar  pressure  distributions  will  be 
exhibited.  Boundary  layer  methods  have  demonstrated  the  capability  to 
calculate  reasonable  boundary  layer  characteristics  by  using  imposed 
pressure  distributions  on  a  flat  plate  calculation  to  model  complex  geo¬ 
metries.  An  attempt  was  made  to  model  the  three-dimensional  duct  flow  of 
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the  Abbiss  test  by  redefining  the  duct  geometry  such  that  the  ratio  of 
total  boundary  layer  displacement  area  to  total  duct  area  was  matched. 

This  approach  assumed  that  a  two-dimensional  duct  with  the  correct  block¬ 
age  area  ratio  distribution  would  produce  the  correct  pressure  distribution 
This  pressure  distribution  would  in  turn  produce  a  correct  boundary  layer 
displacement  thickness  directly  in  the  two-dimensional  calculations. 
Therefore,  once  the  new  duct  height  was  computed  based  on  the  blockage 
area  ratio,  the  modified  duct  geometry  remained  fixed  and  did  not  change 
with  time.  For  the  3  ft  x  3  ft  duct  of  the  Abbiss  test  and  assuming 
negligible  effects  from  the  corners,  the  modified  duct  was  assumed  to  be 
a  1.5  ft  (0.75  ft  to  center  line  boundary)  two-dimensional  duct. 

The  modified  duct  geometry  (I-series)  was  presented  on  Figure  14  as 
the  dash-dot-dash  line.  The  results  show  that  the  x/6^  =  79.5  profile 
was  over-corrected  and  agreement  for  the  downstream  profiles  was  degraded 
relative  to  the  base  line  and  experimental  results.  A  review  of  Figure 
15  showed  that  the  shock  has  been  pushed  forward  and  was  badly  distorted 
at  the  upper  boundary.  A  global  check  of  the  mass  flow  through  each  sub¬ 
sequent  x-mesh  station  revealed  nearly  8%  apparent  mass  was  lost  (in  the 

calculation),  most  of  it  being  lost  at  the  shock  region.  No  conclusion 
about  the  similarity  correction  could  be  made  until  this  large  mass  dis¬ 
crepancy  was  corrected. 

A  similar  mass  flow  check  of  the  base  data  results  in  an  indicated 
2%  mass  gain  at  the  shock  region.  Figure  16  shows  the  shock  for  the  base 
line  case  was  noticeably  distorted  near  the  upper  boundary.  It  was  also 
obvious  that  the  local  distortion  of  the  Mach  contours  (Fig  16a)  was  not 

indicative  of  a  center  line  boundary.  It  was  observed  that  in  cases 

where  mass  was  lower  at  the  outflow  boundary,  the  shock  structure  was 
forced  forward  of  the  desired  shock  location.  Conversely,  where  the  mass 
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Figure  16.  Mach  Contours  Depicting  Improved  Shock-Center  Line  Structure 
With  Corrected  Outflow  Mass  Rate  (0  &  T  -  Series) 
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flow  was  higher  at  the  outflow  boundary  than  at  the  inflow  boundary,  the 
shock  was  distorted  as  if  the  shock  wave  were  being  pushed  downstream  from 
the  fixed  location  on  the  center  line  boundary.  Several  other  cases  were 
checked  and  the  trend  of  these  two  cases  was  consistently  found.  When  the 
mass  flow  was  low,  the  local  velocity  was  low  and  the  local  pressure  was 
high.  Therefore,  a  loss  of  mass  flow  in  the  duct  resulted  in  the  outflow 
pressure  being  too  high. 

An  attempt  to  correct  the  effective  mass  flow  imbalance  was  made. 

The  outflow  boundary  was  assumed  to  be  sufficiently  far  downstream  that 
the  pu  distribution  can  be  assumed  similar  for  the  last  two  streamwise 
mesh  locations.  The  total  mass  flow  was  computed  at  the  inflow  boundary 
and  at  the  ILM1  mesh  station.  The  ratio  of  these  two  totals  was  used  to 
scale  the  ILM1  pu  values  to  compute  outflow  pu  values  which  forced  a 
global  mass  flow  balance  at  the  inflow  and  outflow  boundaries.  This 
condition  replaced  the  base  line  (u)x  =  0  condition  at  the  outflow 
boundary  in  computing  the  T-series  results.  Figure  16b  shows  that  the 
mass  flow  correction  at  the  outflow  boundary  has  almost  totally  eliminated 
the  shock  distortion  at  the  center  line  boundary.  Examination  of  the 
dashed  line  results  (T-series)  in  Figure  14  shows  small  improvements  in 
agreement  with  experimental  results  in  the  near  shock  region  and  an 
almost  equal  degrading  of  agreement  for  the  downstream  profiles.  A  com¬ 
parison  of  the  convergence  plot  for  the  T-series  (Fig  17a)  with  base  line 
plot  (Fig  5,  page  20)  shows  a  three-to-fourfold  increase  in  time  required 
to  reach  a  converged  solution. 

The  last  computation  to  be  performed  was  an  attempt  to  apply  the 
outflow  mass  correction  to  the  modified  duct  three-dimensional  case.  These 
calculations  were  denoted  as  the  U-series.  Figure  17b  shows  that  the  con¬ 
vergence  of  the  U-series  was  far  from  complete  and  Figure  18  reveals  that 
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Figure  18.  Mach  Contours  Depicting  Shock  Movement  with  Attempted 
3-D  and  Mass  Flow  Rate  Corrections  (U-Series) 
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the  shock  was  oscillating  in  the  duct.  However,  Figure  19  indicates 
that  the  flow  was  oscillating  around  the  experimental  data,  This  gives 
encouragement  that  a  three-dimensional  correction  of  this  type  is  feasible. 
The  convergence  of  the  U-series  is  too  slow  to  justify  continued  computa¬ 
tions  without  attempting  to  Improve  the  convergence  rate.  Approximately 
8  hours  of  computer  processor  time  is  represented  on  Figure  17b.  A  very 
rough  damping  ratio  convergence  prediction  calculation  indicated  that  at 
least  another  30  hours  would  be  required  to  converge  the  solution  at  the 
present  rate. 

Conclusions 

The  results  are  summarized  by  the  following  conclusions: 

1.  The  two-dimensional  MacCormack  code  yields  results  that  are  in 
qualitative  and  quantitative  (approximately  10%)  agreement  with  the 
experimental  results  of  Abbiss .  This  agreement  has  been  demonstrated 
for  both  strong  and  weak  transonic  interactions. 

2.  The  experimental  data  exhibits  measurable  three-dimensional 
effects  which  if  accounted  for  in  the  numerical  model  or  elimiated  from 
the  experimental  approach  would  bring  the  experimental  and  numerical 
results  into  even  better  agreement. 

3.  Strong  and  weak  transonic  normal  shock  turbulent  boundary  layer 
interactions  can  be  computed  without  any  modification  to  boundary  con¬ 
dition  treatment. 

4.  The  subsonic  outflow  boundary  treatment  is  very  important.  A 
mathematically  consistent  formulation  to  the  outflow  can  be  physically 
incompatible  with  the  computational  algorithm  and  not  only  affect  the 
solution  by  actually  prevent  convergence. 
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.  5.  The  use  of  a  two-dimensional  code  to  provide  useful  engineering 
results  for  a  three-dimensional  duct  flows  is  both  feasible  and  promising. 

6.  A  more  efficient  method  to  reduce  computing  time  is  very  desirable, 
especially  for  cases  where  three-dimensional  effects  are  required  in  the 
simulation. 


IV.  Recommendations  and  Suggestions 


1.  Further  work  is  needed  relative  to  boundary  cpnditipns  and  their 
numerical  implementation,  especially  for  subsonic  boundaries. 

2.  Effort  should  be  continued  to  determine  a  method  to  implement  the 
shock  wave,  maintain  a  constant  global  mass  flows,  and  insure  a  reasonable 
outflow  pressure  for  this  highly  coupled  problem. 

3.  The  computation  of  the  U-series  (three-dimensional  correction) 
should  be  continued  with  an  improved  convergence  sheme  if  possible. 

4.  More  of  the  experimental  results  should  be  obtained  from  East 
and  Abbiss  and  compared  to  the  numerical  results  to  provide  a  more  com¬ 
plete  analysis  base. 

5.  A  single  case  should  be  computed  with  an  increased  number  of  mesh 
points  (refined  mesh  structure)  to  provide  a  better  check  on  resolution 
effects. 

6.  A  two-dimensional  supersonic  duct  computation  would  provide  a 
better  inflow  boundary  condition  for  the  simulation  than  the  boundary  layer 
code  used  to  date. 

7.  Application  of  other  finite  difference  codes  or  a  finite  element 
technique,  would  provide  another  comparison  with  the  present  methods. 

8.  A  careful  analysis  of  existing  experimental  data  should  be  con¬ 
ducted  to  determine  if  there  is  truly  a  post-shock  supersonic  zone.  A 

new  experimental  effort  may  well  be  required  to  resolve  this  basic  question. 
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Appendix  A 


The  finite  difference  code  used  for  this  study  was  based  on  the 
method  of  MacCormack  (Ref  14).  The  code  is  a  hybrid  implicit-explicit 
scheme  to  solve  the  Navier-Stokes  equations  in  two  dimensions.  The  pro¬ 
gram  is  similar  to  that  employed  by  Shea  (Ref  28) }  but  incorporates 
different  boundary  conditions.  The  procedure  will  be  discussed  in  terms 
of  grid  generation,  numerical  algorithms,  and  turbulence  model.  Boundary 
conditions  are  discussed  separately  in  Sections  II  and  III  with  support¬ 
ing  equation  development  in  Appendix  B. 

The  computational  domain  was  divided  into  a  64  x  32  rectangular 
mesh  system.  The  streamwise  x-direction  was  divided  into  an  equally 
spaced  Ax  mesh  by  64  grid  points.  The  spacing  was  such  that  the  ratio 
of  Ax  to  inflow  boundary  layer  thickness  is  approximately  0.37.  Any 
deviation  from  this  spacing  was  specifically  noted.  The  normal  shock 
wave  was  introduced  into  the  flow  field  between  the  streamwise  grid 
points  16  and  17  where  the  Rankine-Hugonist  conditions  was  applied.  A 
description  of  the  shock  implementation  was  given  in  Section  III. 

The  y-direction  was  divided  into  two  regions.  Points  adjacent  to 
the  wall  were  exponentially  spaced  up  to  approximately  1.2  times  the  inflow 
boundary  layer  thickness.  In  the  outer  region,  the  points  were  equally 
spaced.  Most  of  the  calculations  were  accomplished  with  32  grid  points 
in  the  y-direction  of  which  14  were  in  the  fine  mesh  region.  The  exponen¬ 
tial  spacing  was  generated  by 
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where 


n  (JLFM  -  1) 


So,  ck  was  a  constant  and  n  is  a  nondimens tonal  value  of  y  to  position  the 
J  =  1  and  J  t=  2  grid  points  at  ±0.5  Ay  for  the  reflexion  plane  implementation 
of  the  wall  boundary  conditions.  JLFM  was  the  last  grid  point  in  the  fine 
mesh.  This  relationship  provided  a  mesh  in  which  each  Ay  was  approximately 
0.6,  the  height  of  the  Ay  above  It.  The  value  of  ck  was  adjusted  to  pro¬ 
vide  a  smooth  transition  between  the  fine  mesh  height  and  the  constant  mesh 
si2e  in  the  outer  region.  The  criteria  used  was  that  Ay  in  the  interface 
mesh  (between  J  =  14  and  J  =  15)  was  approximately  the  average  of  the  spac¬ 
ing  just  above  and  below.  This  resulted  in  the  interface  size  of  approxi¬ 
mately  0.8  Ay  of  the  outer  region.  See  Figure  20  for  the  resulting  mesh 
structure. 

The  exponentially  spaced  mesh  adjacent  to  the  wall  allowed  better 
computational  resolution  in  the  viscous  dominated  region  where  high  gra¬ 
dients  occur  along  with  separation  and  reattachment  when  present.  Further, 
the  fine  mesh  was  spaced  so  that  the  first  two  points  at  the  wall  were  kept 
in  the  linear  viscous  subregion  of  the  boundary  layer.  Appendix  B  describes 
the  wall  boundary  implemenation.  The  small  y  increments  close  to  the  wall 
require  consideration  in  computational  algorithm  development. 

The  two-dimensional  time-dependent  compressible  Navier-Stokes  equa¬ 
tions,  A-2,  were  expressed  in  the  flux  vector  form. 

(U)t  *  (F)x  +  (G)y  *  0  A-2 

where 
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Y/DEL(1) 


X/DEL(1) 

a.  Total  Computational  Domain 


X/DKL(1) 

b.  Detailed  Mesh  Structure  Near  the  Wall 


Figure  20,  Computational  Mesh  Structure 


56 


■f 

i 


i 

) 

k 


r  H 


rpu 

2 

pu  +  Ox 

pUV  +  T 
r  xy 

L(e  +  o)u  +  t  V  -  k 
t  x  yx 


G 


+  a 

> 

+  T 


yx 


+  a  )v  +  t 

y  xy 


u 


a  =  p  -  A[(u)  +  (v)  ]  -  zp  (u) 

x  x  y  x 

oy  =  P  -  M(u)x  +  (v)yl  "  ZP  (v)y 

and 

r  =  t  =  -  p[(u)  +  (v)  ] 

xy  yx  y  x 


The  algorithm  used  to  solve  A-2  is  quite  different  for  the  coarse 
and  fine  mesh  regions.  The  fine  mesh  method  was  more  complex  but  provided 
significantly  improved  computing  efficiency  for  the  extremely  small  y 
increments  close  to  the  wall.  For  the  coarse  mesh  the  dependent  variables 
p,  pu,  pv,  and  et  were  estimated  at  a  given  time  t  =  nAt  at  each  grid 
point  i,  j.  The  solution  was  then  computed  at  time  t  =  (n  +  1)  At  by  i 


,n  +  1 

i, 


-  L(At)  uj  . 

■*-»  J 


A- 3 


where  L  was  a  sequence  of  operators  which  explicitly  advanced  the  solution 
At  in  time.  A-3  was  repeated  until  the  change  in  U  between  time  n  and  n  +  1 
was  small  enough  to  assume  the  solution  is  satisfactorily  converged  (See 
Section  III  for  convergence  criteria)  to  a  steady  state  solution.  The 
advance  was  time-split  into  two  primary  operators. 
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L  advances 
x 


(U)t  +  (F)x  "  0 


and  L  advances 

y 


(U)  +  (G)  =0 

y  y 


A-4 


A-5 


The  one  dimensional  difference  operators  L  and  L  were  cast  into  a  sym- 

x  y 


metric  sequence 

n  +  1  .»■  a#-  n 

U,  ,  =  L  (%  L  (At)  L  (~ )  U,  . 


It  > 

.  J  x 


x'  2  i,  j 


A-6 


Thus,  the  operator  was  called  twice  to  advance  the  solution  one  half 

At  each  time  according  to  A-4  with  U  and  F  as  defined  in  A-l.  Between 

the  two  L  steps  of  A-4  L  was  advanced  At  according  to  A-5.  The  combined 
x  y 

effect  of  these  three  sequencial  split  operations  contributed  to  the 
advance  of  A-l  for  a  time  At  (Ref  7). 

The  time  At  was  determined  according  to  the  numerical  stability 
analysis  criteria  of  the  explicit  operators.  The  At  for  each  step  was 
proportional  to  the  x  and  y  mesh  spacing.  In  the  exponentially  spaced 
region,  the  spacing  in  Ay  become  very  small  and  the  allowable  At  was 
proportionally  small.  Therefore,  to  improve  efficiency  a  more  complex 
algorithm  was  used  in  this  region. 

The  fine  mesh  algorithm  further  subdivided  L  into  its  L  inviscid 


H 


(hyperbolic)  terms  and  L  viscous  (parabolic)  terms. 

yp 

L  advances 


'H 


(U)  +  (G  )  *  0 

y  h  v 


A-7 


where 


'H  * 
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by  an  explicit  approach  using  characteristic  relations  to  predict  con¬ 
vection  and  pressure  fields  and  was  stable  proyided 


At  <  i^i 
~  |V  I 


This  resulted  in  a  much  larger  allowable  time  step  because  as  Ay 
becomes  very  small  close  to  the  wall  v  is  also  approaching  zero.  The 
result  was  stable  time  steps  which  were  50  to  100  times  larger  than  other¬ 
wise  possible  (Ref  14:7). 

L  advances 


where 


<U)  +  (G  )  =  0 

t  p  y 


G  =  G  -  G 
P  H 


A-8 


by  an  implicit  procedure  of  either  Crank-Nicolson  or  Laasonen  type.  The 

2 

second-order  accurate  Crank-Nicolson  method  was  used  when  pAt/dy  was  small 
enough  to  preclude  erratic  behavior.  When  this  quantity  was  large  enough 
to  cause  numerical  instability,  the  first  order  Laasonen  method  was  used. 
The  L  and  L  operators  were  incorporated  in  the  fine  mesh  as 


U 


P 

n  +  1 

i*  j 


'H 


l  L  (4^)  L  (I1)  L  (— )  L  .At.  /At™  “ 

y„  2ra  y  2m  x  m  y  (-r— )  L  (■=—)]  U. 

H  Jp  '2n/  yHv2m/J  i,  j 


where  m  was  the  number  of  fine  mesh  time  step  required  to  equal  one  course 
mesh  time  step.  Thus,  the  fine  mesh  was  advanced  by  five  sequential 
operators  m  times  every  time  the  coarse  mesh  was  advanced  by  time  At. 

All  of  the  operators  were  conservative  and  At  was  computed  for  each 
step  to  insure  numerical  stability.  Ref  14  provided  a  much  more  complete 
description  of  the  numerical  operators  used  for  this  effort. 
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The  turbulence  model  used  was  the  algebraic  eddy  viscosity  model  used 
by  Shea  (Ref  28).  It  was  a  modified  Cebeci-Smith-Mosinakis  model  in  the 
inner  layer  and  a  Maise-McDonald  model  in  the  law  of  the  wake  region. 
References  6,  30,  and  31  indicated  that  this  model  was  probably  most  appro¬ 
priate  and  the  incorporation  of  more  complex  schemes  would  not  result  in 
improved  solutions. 
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Appendix  B 

Boundary  Condition  Supporting  Equations 

Wall  Boundary  Finite  Difference  Equations 

The  reflection  plane  implementation  at  the  wall  of. 


v  =  u  =  (T)  “  0  B-l 

y 

with  J  *  1  at  -0.5  Ay(l)  and  J  =  2  at  +0.5  A y(l)  from  the  wall  was: 

u(l,  1)  =  — u(I,  2)  .  B-2 

v(I,  1)  =  -v(l,  2)  B-3 

El (I,  1)  =  EI(I,  2)  B-4 

where  El  denotes  specific  internal  energy. 

The  compatible  condition  for  the  normal  (y)  momentum  equation  at  the 


wall  Eqn  11-3  was  implemented  (after  Ref  26)  to  provide  pressure  at  the 
wall  by: 

Pi,  1  =  Pi,  2 

-  Aytpi+l(uI+1.3  ~  ui+l . 1^  "  vi-l(uj-1.3  ~  ui-l.l^ 

2Ax(y3  -  yx) 

+  XiI(Ui+1.2  ~  “i-1.2*  "  (ui+l.l  ~  Ui-l.l)] 

2  Ax  Ay  g_c 

(ti.3  -  Va*  <vi.2  - 

>3  -  >2 _ _ 

o*5(y3  "  y{i 
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+  <xi  +  2vi> 


W 


Center  Line  Boundary 


The  steady  2-D  governing  equations  assuming  the  flow 

to  be  inviscid 

and  adiabatic  at  the  center  line  were: 

Continuity 

(pu)x  +  ^Pv)y  “  0 

B-6 

x-momentum 

(pu2)x  +  (puv)y  +  (p)x  «*  0 

B-7 

y-momentum 

(puv)x  +  (pv2)y  +  (p)y=  0 

B-8 

energy 

(puet)x  +  (pvet)y  =  0 

B-9 

Applying  continuity  (B-6)  to  B-7,  8,  and  9  yields: 

pu(u)x  +  pv(u)y  +  (p)x  =  0 

B-10 

Pu(v)  +  pv(v)  +  (p)  =  0 

a  y  y 

B-ll 

pu(e. )  +  pv(e „)  ®  0 

t  x  K  t  y 

B-12 

Introducing  the  v=  0  condition  at  the  center  into  B-6,  10, 

,  11  and  12 

yield: 

(pu)x  +  <v)y  =  0 

B-13 

pu(u)x  +  (p)x  =  0 

B-14 

(P)y  »  0 

B-15 

‘V*  ■ 0 

B-16 

Therefore,  (et>x  ■  0  can  be  added  to  the  list  of  compatible  center  line 
boundary  conditions. 
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20.  ABSTRACT  (Continue  on  reeetee  tide  II  neceeeery  and  Identity  by  block  number) 

The  hybrid  finite  difference  code  developed  by  MacCormack  was  applied  to  the 
investigation  of  transonic  normal-shock  turbulent  boundary  layer  interactions. 
The  computations  were  performed  for  the  half  plane  of  a  symmetric  two  dimen¬ 
sional  duct  by  establishing  a  symmetry  boundary  condition  at  the  upper  bound¬ 
ary.  Both  first  and  second  order  center  line  boundary  conditions  were  imposed 
with  no  measurable  difference  observed.  A  two-point  linear  extrapolation  of 
the.prlmative  variable  was  unsuccessfully  attempted  at  the  subsonic  outflow 
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boundary,  but  a  simple  zero  gradient  condition  gave  satisfactory  results  at 
four  different  outflow  boundary  positions  relative  to  the  shock  wave.  Numeri¬ 
cal  results  (M  -  1.51,  1.40  and  1.30  Re  -  3X10®  per  ft)  were  compared  with 
the  experimental  data  reported  by  Abblss  and  East.  Even  though  the  data 
exhibit  three-dimensional  effects,  the  two-dimensional  computations  show 
agreement  within  approximately  10%.  The  difference  observed  In  the  numerical- 
experimental  comparisons  were  all  consistent  with  expected  three-dimensional 
trends.  Although  not  conclusive,  the  potential  of  adding  simple  three-dlmen- 
sionalcorrections  to  the  two-dimensional  code  shows  promise  for  Improving  the 
experimental-numerical  agreement . 


